{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import linecache\n",
    "import datetime\n",
    "import sys\n",
    "import numpy as np\n",
    "import random\n",
    "from scipy.sparse import *\n",
    "from scipy import *\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import os.path\n",
    "import matplotlib\n",
    "import copy\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from tempfile import TemporaryFile\n",
    "import pickle\n",
    "def flushPrint(variable):\n",
    "    sys.stdout.write('\\r')\n",
    "    sys.stdout.write('%s' % variable)\n",
    "    sys.stdout.flush()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 读取2016年的流量数据并与2020最新部分比例数据合并"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pi</th>\n",
       "      <th>pj</th>\n",
       "      <th>loni</th>\n",
       "      <th>lati</th>\n",
       "      <th>cityi.id</th>\n",
       "      <th>provi.id</th>\n",
       "      <th>lonj</th>\n",
       "      <th>latj</th>\n",
       "      <th>cityj.id</th>\n",
       "      <th>provj.id</th>\n",
       "      <th>flowij</th>\n",
       "      <th>distij</th>\n",
       "      <th>sameprov</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1220637</td>\n",
       "      <td>1220637</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>1082923</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1220637</td>\n",
       "      <td>488412</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>117.344841</td>\n",
       "      <td>39.284732</td>\n",
       "      <td>120000</td>\n",
       "      <td>120000</td>\n",
       "      <td>4642</td>\n",
       "      <td>128.019884</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1220637</td>\n",
       "      <td>371707</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>114.439494</td>\n",
       "      <td>38.130579</td>\n",
       "      <td>130100</td>\n",
       "      <td>130000</td>\n",
       "      <td>2962</td>\n",
       "      <td>284.863566</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1220637</td>\n",
       "      <td>188880</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>118.335459</td>\n",
       "      <td>39.713888</td>\n",
       "      <td>130200</td>\n",
       "      <td>130000</td>\n",
       "      <td>1632</td>\n",
       "      <td>172.096010</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1220637</td>\n",
       "      <td>82307</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>119.187239</td>\n",
       "      <td>40.085314</td>\n",
       "      <td>130300</td>\n",
       "      <td>130000</td>\n",
       "      <td>1021</td>\n",
       "      <td>236.130974</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1220637</td>\n",
       "      <td>214644</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>114.542925</td>\n",
       "      <td>36.552456</td>\n",
       "      <td>130400</td>\n",
       "      <td>130000</td>\n",
       "      <td>3332</td>\n",
       "      <td>435.600363</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1220637</td>\n",
       "      <td>138234</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>114.816896</td>\n",
       "      <td>37.212320</td>\n",
       "      <td>130500</td>\n",
       "      <td>130000</td>\n",
       "      <td>1873</td>\n",
       "      <td>358.424618</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1220637</td>\n",
       "      <td>293124</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>115.171130</td>\n",
       "      <td>39.021666</td>\n",
       "      <td>130600</td>\n",
       "      <td>130000</td>\n",
       "      <td>6091</td>\n",
       "      <td>167.515139</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1220637</td>\n",
       "      <td>87130</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>115.031804</td>\n",
       "      <td>40.864962</td>\n",
       "      <td>130700</td>\n",
       "      <td>130000</td>\n",
       "      <td>2891</td>\n",
       "      <td>139.016370</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1220637</td>\n",
       "      <td>61064</td>\n",
       "      <td>116.412574</td>\n",
       "      <td>40.185609</td>\n",
       "      <td>110000</td>\n",
       "      <td>110000</td>\n",
       "      <td>117.547029</td>\n",
       "      <td>41.347208</td>\n",
       "      <td>130800</td>\n",
       "      <td>130000</td>\n",
       "      <td>1866</td>\n",
       "      <td>160.655042</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        pi       pj        loni       lati  cityi.id  provi.id        lonj  \\\n",
       "0  1220637  1220637  116.412574  40.185609    110000    110000  116.412574   \n",
       "1  1220637   488412  116.412574  40.185609    110000    110000  117.344841   \n",
       "2  1220637   371707  116.412574  40.185609    110000    110000  114.439494   \n",
       "3  1220637   188880  116.412574  40.185609    110000    110000  118.335459   \n",
       "4  1220637    82307  116.412574  40.185609    110000    110000  119.187239   \n",
       "5  1220637   214644  116.412574  40.185609    110000    110000  114.542925   \n",
       "6  1220637   138234  116.412574  40.185609    110000    110000  114.816896   \n",
       "7  1220637   293124  116.412574  40.185609    110000    110000  115.171130   \n",
       "8  1220637    87130  116.412574  40.185609    110000    110000  115.031804   \n",
       "9  1220637    61064  116.412574  40.185609    110000    110000  117.547029   \n",
       "\n",
       "        latj  cityj.id  provj.id   flowij      distij  sameprov  \n",
       "0  40.185609    110000    110000  1082923    0.000000         1  \n",
       "1  39.284732    120000    120000     4642  128.019884         0  \n",
       "2  38.130579    130100    130000     2962  284.863566         0  \n",
       "3  39.713888    130200    130000     1632  172.096010         0  \n",
       "4  40.085314    130300    130000     1021  236.130974         0  \n",
       "5  36.552456    130400    130000     3332  435.600363         0  \n",
       "6  37.212320    130500    130000     1873  358.424618         0  \n",
       "7  39.021666    130600    130000     6091  167.515139         0  \n",
       "8  40.864962    130700    130000     2891  139.016370         0  \n",
       "9  41.347208    130800    130000     1866  160.655042         0  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#读取工作簿和工作簿中的工作表\n",
    "df=pd.read_csv('city_flow_v1.csv')\n",
    "df.head(10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.995342957403680654"
     ]
    }
   ],
   "source": [
    "cityi = list(df['cityi.id'])\n",
    "cityj = list(df['cityj.id'])\n",
    "all_cities = set(cityi+cityj)\n",
    "cities_temp = {}\n",
    "flows = {}\n",
    "for n,i in enumerate(df.iterrows()):\n",
    "    if n % 1000 == 0:\n",
    "        flushPrint(n/len(df))\n",
    "    cityi=long(i[1]['cityi.id'])\n",
    "    cityj=long(i[1]['cityj.id'])\n",
    "    provi=long(i[1]['provi.id'])\n",
    "    provj=long(i[1]['provj.id'])\n",
    "    cities_temp[cityi] = {'id':cityi,'lon':i[1]['loni'],'lat':i[1]['lati'],'provid':provi,'pop':i[1]['pi']}\n",
    "    cities_temp[cityj] = {'id':cityj,'lon':i[1]['lonj'],'lat':i[1]['latj'],'provid':provj,'pop':i[1]['pj']}\n",
    "    value = flows.get((cityi,cityj),0)\n",
    "    flows[(cityi,cityj)] = value + i[1]['flowij']\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#存储到流量矩阵matrix\n",
    "nodes_temp={} #cityid: row index mapping\n",
    "matrix = np.zeros([len(cities_temp), len(cities_temp)])\n",
    "self_flux = np.zeros(len(cities_temp))\n",
    "pij1 = np.zeros([len(cities_temp), len(cities_temp)])\n",
    "for key,value in flows.items():\n",
    "    id1 = nodes_temp.get(key[0],-1)\n",
    "    if id1<0:\n",
    "        id1 = len(nodes_temp)\n",
    "        nodes_temp[key[0]] = len(nodes_temp)\n",
    "    id2 = nodes_temp.get(key[1],-1)\n",
    "    if id2<0:\n",
    "        id2 = len(nodes_temp)\n",
    "        nodes_temp[key[1]] = len(nodes_temp)\n",
    "    matrix[id1, id2] = value\n",
    "for i in range(matrix.shape[0]):\n",
    "    self_flux[i] = matrix[i, i]\n",
    "    matrix[i, i] = 0\n",
    "    if np.sum(matrix[i,:])>0:\n",
    "        pij1[i,:]=matrix[i,:]/np.sum(matrix[i,:])\n",
    "#根据我国人口总数，核算每个城市的人口数以及总流动人口占比\n",
    "all_pop = np.sum(np.array([v['pop'] for k,v in cities_temp.items()]))\n",
    "china_pop = 13*1e8\n",
    "frac = china_pop / all_pop\n",
    "for k,v in cities_temp.items():\n",
    "    cities_temp[k]['pop'] = frac * v['pop']\n",
    "flowing_ratio=np.sum(matrix) / all_pop\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 读取每个城市id与name的映射"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>COUNTY_ID</th>\n",
       "      <th>COUNTY</th>\n",
       "      <th>CITY_ID</th>\n",
       "      <th>CITY</th>\n",
       "      <th>PROV_ID</th>\n",
       "      <th>PROV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>110101</td>\n",
       "      <td>东城区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>110102</td>\n",
       "      <td>西城区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>110105</td>\n",
       "      <td>朝阳区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>110106</td>\n",
       "      <td>丰台区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>110107</td>\n",
       "      <td>石景山区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>110108</td>\n",
       "      <td>海淀区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>110109</td>\n",
       "      <td>门头沟区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>110111</td>\n",
       "      <td>房山区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>110112</td>\n",
       "      <td>通州区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>110113</td>\n",
       "      <td>顺义区</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "      <td>110000</td>\n",
       "      <td>北京市</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   COUNTY_ID COUNTY  CITY_ID CITY  PROV_ID PROV\n",
       "0     110101    东城区   110000  北京市   110000  北京市\n",
       "1     110102    西城区   110000  北京市   110000  北京市\n",
       "2     110105    朝阳区   110000  北京市   110000  北京市\n",
       "3     110106    丰台区   110000  北京市   110000  北京市\n",
       "4     110107   石景山区   110000  北京市   110000  北京市\n",
       "5     110108    海淀区   110000  北京市   110000  北京市\n",
       "6     110109   门头沟区   110000  北京市   110000  北京市\n",
       "7     110111    房山区   110000  北京市   110000  北京市\n",
       "8     110112    通州区   110000  北京市   110000  北京市\n",
       "9     110113    顺义区   110000  北京市   110000  北京市"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df=pd.read_csv('county_city_province.csv')\n",
    "df.head(10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_to_city={}\n",
    "for i in df.iterrows():\n",
    "    id_to_city[i[1]['CITY_ID']]={'city':i[1]['CITY'],'prov':i[1]['PROV']}\n",
    "city_properties = {}\n",
    "nodes = {}\n",
    "for k,v in cities_temp.items():\n",
    "    new_property = id_to_city.get(k, {'city':str(k), 'prov':str(k)})\n",
    "    v['prov'] = new_property['prov']\n",
    "    city_properties[new_property['city']] = v\n",
    "    row_idx = nodes_temp.get(v['id'],-1)\n",
    "    if row_idx >= 0:\n",
    "        nodes[new_property['city']] = row_idx\n",
    "    else:\n",
    "        print(new_property['city'])\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'北京市': 0,\n",
       " '天津市': 1,\n",
       " '石家庄市': 2,\n",
       " '唐山市': 3,\n",
       " '秦皇岛市': 4,\n",
       " '邯郸市': 5,\n",
       " '邢台市': 6,\n",
       " '保定市': 7,\n",
       " '张家口市': 8,\n",
       " '承德市': 9,\n",
       " '沧州市': 10,\n",
       " '廊坊市': 11,\n",
       " '衡水市': 12,\n",
       " '太原市': 13,\n",
       " '大同市': 14,\n",
       " '阳泉市': 15,\n",
       " '长治市': 16,\n",
       " '晋城市': 17,\n",
       " '朔州市': 18,\n",
       " '晋中市': 19,\n",
       " '运城市': 20,\n",
       " '忻州市': 21,\n",
       " '临汾市': 22,\n",
       " '吕梁市': 23,\n",
       " '呼和浩特市': 24,\n",
       " '包头市': 25,\n",
       " '乌海市': 26,\n",
       " '赤峰市': 27,\n",
       " '通辽市': 28,\n",
       " '鄂尔多斯市': 29,\n",
       " '呼伦贝尔市': 30,\n",
       " '巴彦淖尔市': 31,\n",
       " '乌兰察布市': 32,\n",
       " '兴安盟': 33,\n",
       " '锡林郭勒盟': 34,\n",
       " '阿拉善盟': 35,\n",
       " '沈阳市': 36,\n",
       " '大连市': 37,\n",
       " '鞍山市': 38,\n",
       " '抚顺市': 39,\n",
       " '本溪市': 40,\n",
       " '丹东市': 41,\n",
       " '锦州市': 42,\n",
       " '营口市': 43,\n",
       " '阜新市': 44,\n",
       " '辽阳市': 45,\n",
       " '盘锦市': 46,\n",
       " '铁岭市': 47,\n",
       " '朝阳市': 48,\n",
       " '葫芦岛市': 49,\n",
       " '长春市': 50,\n",
       " '吉林市': 51,\n",
       " '四平市': 52,\n",
       " '辽源市': 53,\n",
       " '通化市': 54,\n",
       " '白山市': 55,\n",
       " '松原市': 56,\n",
       " '白城市': 57,\n",
       " '延边朝鲜族自治州': 58,\n",
       " '哈尔滨市': 59,\n",
       " '齐齐哈尔市': 60,\n",
       " '鸡西市': 61,\n",
       " '鹤岗市': 62,\n",
       " '双鸭山市': 63,\n",
       " '大庆市': 64,\n",
       " '伊春市': 65,\n",
       " '佳木斯市': 66,\n",
       " '七台河市': 67,\n",
       " '牡丹江市': 68,\n",
       " '黑河市': 69,\n",
       " '绥化市': 70,\n",
       " '大兴安岭地区': 71,\n",
       " '上海市': 72,\n",
       " '南京市': 73,\n",
       " '无锡市': 74,\n",
       " '徐州市': 75,\n",
       " '常州市': 76,\n",
       " '苏州市': 77,\n",
       " '南通市': 78,\n",
       " '连云港市': 79,\n",
       " '淮安市': 80,\n",
       " '盐城市': 81,\n",
       " '扬州市': 82,\n",
       " '镇江市': 83,\n",
       " '泰州市': 84,\n",
       " '宿迁市': 85,\n",
       " '杭州市': 86,\n",
       " '宁波市': 87,\n",
       " '温州市': 88,\n",
       " '嘉兴市': 89,\n",
       " '湖州市': 90,\n",
       " '绍兴市': 91,\n",
       " '金华市': 92,\n",
       " '衢州市': 93,\n",
       " '舟山市': 94,\n",
       " '台州市': 95,\n",
       " '丽水市': 96,\n",
       " '合肥市': 97,\n",
       " '芜湖市': 98,\n",
       " '蚌埠市': 99,\n",
       " '淮南市': 100,\n",
       " '马鞍山市': 101,\n",
       " '淮北市': 102,\n",
       " '铜陵市': 103,\n",
       " '安庆市': 104,\n",
       " '黄山市': 105,\n",
       " '滁州市': 106,\n",
       " '阜阳市': 107,\n",
       " '341300': 108,\n",
       " '六安市': 109,\n",
       " '亳州市': 110,\n",
       " '池州市': 111,\n",
       " '宣城市': 112,\n",
       " '福州市': 113,\n",
       " '厦门市': 114,\n",
       " '莆田市': 115,\n",
       " '三明市': 116,\n",
       " '泉州市': 117,\n",
       " '漳州市': 118,\n",
       " '南平市': 119,\n",
       " '龙岩市': 120,\n",
       " '宁德市': 121,\n",
       " '南昌市': 122,\n",
       " '景德镇市': 123,\n",
       " '萍乡市': 124,\n",
       " '九江市': 125,\n",
       " '新余市': 126,\n",
       " '鹰潭市': 127,\n",
       " '赣州市': 128,\n",
       " '吉安市': 129,\n",
       " '宜春市': 130,\n",
       " '抚州市': 131,\n",
       " '上饶市': 132,\n",
       " '济南市': 133,\n",
       " '青岛市': 134,\n",
       " '淄博市': 135,\n",
       " '枣庄市': 136,\n",
       " '东营市': 137,\n",
       " '烟台市': 138,\n",
       " '潍坊市': 139,\n",
       " '济宁市': 140,\n",
       " '泰安市': 141,\n",
       " '威海市': 142,\n",
       " '日照市': 143,\n",
       " '莱芜市': 144,\n",
       " '临沂市': 145,\n",
       " '德州市': 146,\n",
       " '聊城市': 147,\n",
       " '滨州市': 148,\n",
       " '菏泽市': 149,\n",
       " '郑州市': 150,\n",
       " '开封市': 151,\n",
       " '洛阳市': 152,\n",
       " '平顶山市': 153,\n",
       " '安阳市': 154,\n",
       " '鹤壁市': 155,\n",
       " '新乡市': 156,\n",
       " '焦作市': 157,\n",
       " '濮阳市': 158,\n",
       " '许昌市': 159,\n",
       " '漯河市': 160,\n",
       " '三门峡市': 161,\n",
       " '南阳市': 162,\n",
       " '商丘市': 163,\n",
       " '信阳市': 164,\n",
       " '周口市': 165,\n",
       " '驻马店市': 166,\n",
       " '济源市': 167,\n",
       " '武汉市': 168,\n",
       " '黄石市': 169,\n",
       " '十堰市': 170,\n",
       " '宜昌市': 171,\n",
       " '襄阳市': 172,\n",
       " '鄂州市': 173,\n",
       " '荆门市': 174,\n",
       " '孝感市': 175,\n",
       " '荆州市': 176,\n",
       " '黄冈市': 177,\n",
       " '咸宁市': 178,\n",
       " '随州市': 179,\n",
       " '恩施土家族苗族自治州': 180,\n",
       " '仙桃市': 181,\n",
       " '潜江市': 182,\n",
       " '天门市': 183,\n",
       " '神农架林区': 184,\n",
       " '长沙市': 185,\n",
       " '株洲市': 186,\n",
       " '湘潭市': 187,\n",
       " '衡阳市': 188,\n",
       " '邵阳市': 189,\n",
       " '岳阳市': 190,\n",
       " '常德市': 191,\n",
       " '张家界市': 192,\n",
       " '益阳市': 193,\n",
       " '郴州市': 194,\n",
       " '永州市': 195,\n",
       " '怀化市': 196,\n",
       " '娄底市': 197,\n",
       " '湘西土家族苗族自治州': 198,\n",
       " '广州市': 199,\n",
       " '韶关市': 200,\n",
       " '深圳市': 201,\n",
       " '珠海市': 202,\n",
       " '汕头市': 203,\n",
       " '佛山市': 204,\n",
       " '江门市': 205,\n",
       " '湛江市': 206,\n",
       " '茂名市': 207,\n",
       " '肇庆市': 208,\n",
       " '惠州市': 209,\n",
       " '梅州市': 210,\n",
       " '汕尾市': 211,\n",
       " '河源市': 212,\n",
       " '阳江市': 213,\n",
       " '清远市': 214,\n",
       " '东莞市': 215,\n",
       " '中山市': 216,\n",
       " '东沙群岛': 217,\n",
       " '潮州市': 218,\n",
       " '揭阳市': 219,\n",
       " '云浮市': 220,\n",
       " '南宁市': 221,\n",
       " '柳州市': 222,\n",
       " '桂林市': 223,\n",
       " '梧州市': 224,\n",
       " '北海市': 225,\n",
       " '防城港市': 226,\n",
       " '钦州市': 227,\n",
       " '贵港市': 228,\n",
       " '玉林市': 229,\n",
       " '百色市': 230,\n",
       " '贺州市': 231,\n",
       " '河池市': 232,\n",
       " '来宾市': 233,\n",
       " '崇左市': 234,\n",
       " '海口市': 235,\n",
       " '三亚市': 236,\n",
       " '三沙市': 237,\n",
       " '儋州市': 238,\n",
       " '五指山市': 239,\n",
       " '琼海市': 240,\n",
       " '文昌市': 241,\n",
       " '万宁市': 242,\n",
       " '东方市': 243,\n",
       " '定安县': 244,\n",
       " '屯昌县': 245,\n",
       " '澄迈县': 246,\n",
       " '临高县': 247,\n",
       " '白沙黎族自治县': 248,\n",
       " '昌江黎族自治县': 249,\n",
       " '乐东黎族自治县': 250,\n",
       " '陵水黎族自治县': 251,\n",
       " '保亭黎族苗族自治县': 252,\n",
       " '琼中黎族苗族自治县': 253,\n",
       " '重庆市': 254,\n",
       " '成都市': 255,\n",
       " '自贡市': 256,\n",
       " '攀枝花市': 257,\n",
       " '泸州市': 258,\n",
       " '德阳市': 259,\n",
       " '绵阳市': 260,\n",
       " '广元市': 261,\n",
       " '遂宁市': 262,\n",
       " '内江市': 263,\n",
       " '乐山市': 264,\n",
       " '南充市': 265,\n",
       " '眉山市': 266,\n",
       " '宜宾市': 267,\n",
       " '广安市': 268,\n",
       " '达州市': 269,\n",
       " '雅安市': 270,\n",
       " '巴中市': 271,\n",
       " '资阳市': 272,\n",
       " '阿坝藏族羌族自治州': 273,\n",
       " '甘孜藏族自治州': 274,\n",
       " '凉山彝族自治州': 275,\n",
       " '贵阳市': 276,\n",
       " '六盘水市': 277,\n",
       " '遵义市': 278,\n",
       " '安顺市': 279,\n",
       " '毕节市': 280,\n",
       " '铜仁市': 281,\n",
       " '黔西南布依族苗族自治州': 282,\n",
       " '黔东南苗族侗族自治州': 283,\n",
       " '黔南布依族苗族自治州': 284,\n",
       " '昆明市': 285,\n",
       " '曲靖市': 286,\n",
       " '玉溪市': 287,\n",
       " '保山市': 288,\n",
       " '昭通市': 289,\n",
       " '丽江市': 290,\n",
       " '普洱市': 291,\n",
       " '临沧市': 292,\n",
       " '楚雄彝族自治州': 293,\n",
       " '红河哈尼族彝族自治州': 294,\n",
       " '文山壮族苗族自治州': 295,\n",
       " '西双版纳傣族自治州': 296,\n",
       " '大理白族自治州': 297,\n",
       " '德宏傣族景颇族自治州': 298,\n",
       " '怒江傈僳族自治州': 299,\n",
       " '迪庆藏族自治州': 300,\n",
       " '拉萨市': 301,\n",
       " '日喀则市': 302,\n",
       " '昌都市': 303,\n",
       " '林芝市': 304,\n",
       " '山南市': 305,\n",
       " '那曲地区': 306,\n",
       " '阿里地区': 307,\n",
       " '西安市': 308,\n",
       " '铜川市': 309,\n",
       " '宝鸡市': 310,\n",
       " '咸阳市': 311,\n",
       " '渭南市': 312,\n",
       " '延安市': 313,\n",
       " '汉中市': 314,\n",
       " '榆林市': 315,\n",
       " '安康市': 316,\n",
       " '商洛市': 317,\n",
       " '兰州市': 318,\n",
       " '嘉峪关市': 319,\n",
       " '金昌市': 320,\n",
       " '白银市': 321,\n",
       " '天水市': 322,\n",
       " '武威市': 323,\n",
       " '张掖市': 324,\n",
       " '平凉市': 325,\n",
       " '酒泉市': 326,\n",
       " '庆阳市': 327,\n",
       " '定西市': 328,\n",
       " '陇南市': 329,\n",
       " '临夏回族自治州': 330,\n",
       " '甘南藏族自治州': 331,\n",
       " '西宁市': 332,\n",
       " '海东市': 333,\n",
       " '海北藏族自治州': 334,\n",
       " '黄南藏族自治州': 335,\n",
       " '海南藏族自治州': 336,\n",
       " '果洛藏族自治州': 337,\n",
       " '玉树藏族自治州': 338,\n",
       " '海西蒙古族藏族自治州': 339,\n",
       " '银川市': 340,\n",
       " '石嘴山市': 341,\n",
       " '吴忠市': 342,\n",
       " '固原市': 343,\n",
       " '中卫市': 344,\n",
       " '乌鲁木齐市': 345,\n",
       " '克拉玛依市': 346,\n",
       " '吐鲁番市': 347,\n",
       " '哈密市': 348,\n",
       " '昌吉回族自治州': 349,\n",
       " '博尔塔拉蒙古族自治州': 350,\n",
       " '巴音郭楞蒙古自治州': 351,\n",
       " '阿克苏地区': 352,\n",
       " '克孜勒苏柯尔克孜自治州': 353,\n",
       " '喀什地区': 354,\n",
       " '和田地区': 355,\n",
       " '伊犁哈萨克自治州': 356,\n",
       " '塔城地区': 357,\n",
       " '阿勒泰地区': 358,\n",
       " '石河子市': 359,\n",
       " '阿拉尔市': 360,\n",
       " '图木舒克市': 361,\n",
       " '五家渠市': 362,\n",
       " '北屯市': 363,\n",
       " '铁门关市': 364,\n",
       " '双河市': 365,\n",
       " '可克达拉市': 366,\n",
       " '昆玉市': 367,\n",
       " '台湾省': 368,\n",
       " '香港特别行政区': 369,\n",
       " '澳门特别行政区': 370}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 读取2020年的百度部分流量比例数据 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Cities</th>\n",
       "      <th>北京市</th>\n",
       "      <th>天津市</th>\n",
       "      <th>石家庄市</th>\n",
       "      <th>唐山市</th>\n",
       "      <th>秦皇岛市</th>\n",
       "      <th>邯郸市</th>\n",
       "      <th>邢台市</th>\n",
       "      <th>保定市</th>\n",
       "      <th>张家口市</th>\n",
       "      <th>...</th>\n",
       "      <th>和田地区</th>\n",
       "      <th>伊犁哈萨克自治州</th>\n",
       "      <th>塔城地区</th>\n",
       "      <th>阿勒泰地区</th>\n",
       "      <th>石河子市</th>\n",
       "      <th>阿拉尔市</th>\n",
       "      <th>五家渠市</th>\n",
       "      <th>香港特别行政区</th>\n",
       "      <th>澳门特别行政区</th>\n",
       "      <th>Unnamed: 352</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>北京市</td>\n",
       "      <td>inf</td>\n",
       "      <td>11.855277</td>\n",
       "      <td>5.725960</td>\n",
       "      <td>12.749903</td>\n",
       "      <td>9.573570</td>\n",
       "      <td>6.598759</td>\n",
       "      <td>4.861948</td>\n",
       "      <td>27.363624</td>\n",
       "      <td>37.854040</td>\n",
       "      <td>...</td>\n",
       "      <td>0.729644</td>\n",
       "      <td>0.424284</td>\n",
       "      <td>0.146950</td>\n",
       "      <td>0.244898</td>\n",
       "      <td>0.284529</td>\n",
       "      <td>0.192485</td>\n",
       "      <td>0.109344</td>\n",
       "      <td>1.716803</td>\n",
       "      <td>1.431816</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>天津市</td>\n",
       "      <td>4.458202</td>\n",
       "      <td>inf</td>\n",
       "      <td>1.806775</td>\n",
       "      <td>20.849449</td>\n",
       "      <td>5.471267</td>\n",
       "      <td>5.470594</td>\n",
       "      <td>1.854504</td>\n",
       "      <td>3.523714</td>\n",
       "      <td>3.469709</td>\n",
       "      <td>...</td>\n",
       "      <td>0.527886</td>\n",
       "      <td>0.085204</td>\n",
       "      <td>0.012593</td>\n",
       "      <td>0.069344</td>\n",
       "      <td>0.005139</td>\n",
       "      <td>0.027000</td>\n",
       "      <td>0.026867</td>\n",
       "      <td>0.228988</td>\n",
       "      <td>0.325382</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>石家庄市</td>\n",
       "      <td>2.230990</td>\n",
       "      <td>1.815094</td>\n",
       "      <td>inf</td>\n",
       "      <td>3.832100</td>\n",
       "      <td>3.980366</td>\n",
       "      <td>11.405949</td>\n",
       "      <td>30.187961</td>\n",
       "      <td>15.307815</td>\n",
       "      <td>5.608304</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.081481</td>\n",
       "      <td>0.003760</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.075537</td>\n",
       "      <td>0.006981</td>\n",
       "      <td>0.044484</td>\n",
       "      <td>0.012442</td>\n",
       "      <td>0.007495</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>唐山市</td>\n",
       "      <td>1.755526</td>\n",
       "      <td>6.625320</td>\n",
       "      <td>2.427825</td>\n",
       "      <td>inf</td>\n",
       "      <td>23.561492</td>\n",
       "      <td>1.150411</td>\n",
       "      <td>0.838385</td>\n",
       "      <td>1.560130</td>\n",
       "      <td>2.586040</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.012455</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.074082</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>秦皇岛市</td>\n",
       "      <td>0.743911</td>\n",
       "      <td>1.104723</td>\n",
       "      <td>0.999676</td>\n",
       "      <td>12.371138</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.347803</td>\n",
       "      <td>0.232399</td>\n",
       "      <td>0.585399</td>\n",
       "      <td>0.724584</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>邯郸市</td>\n",
       "      <td>3.058310</td>\n",
       "      <td>3.685172</td>\n",
       "      <td>9.266465</td>\n",
       "      <td>2.033850</td>\n",
       "      <td>1.596207</td>\n",
       "      <td>inf</td>\n",
       "      <td>18.648848</td>\n",
       "      <td>2.832027</td>\n",
       "      <td>1.373215</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.010470</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.040079</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>邢台市</td>\n",
       "      <td>1.784559</td>\n",
       "      <td>1.592507</td>\n",
       "      <td>17.716423</td>\n",
       "      <td>1.595459</td>\n",
       "      <td>1.481331</td>\n",
       "      <td>17.228512</td>\n",
       "      <td>inf</td>\n",
       "      <td>2.246571</td>\n",
       "      <td>1.065858</td>\n",
       "      <td>...</td>\n",
       "      <td>0.009354</td>\n",
       "      <td>0.031075</td>\n",
       "      <td>0.004738</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.034814</td>\n",
       "      <td>0.009440</td>\n",
       "      <td>0.009669</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>保定市</td>\n",
       "      <td>7.296308</td>\n",
       "      <td>3.120498</td>\n",
       "      <td>14.208043</td>\n",
       "      <td>3.131967</td>\n",
       "      <td>2.768861</td>\n",
       "      <td>2.910327</td>\n",
       "      <td>2.589855</td>\n",
       "      <td>inf</td>\n",
       "      <td>5.916331</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.007716</td>\n",
       "      <td>0.005211</td>\n",
       "      <td>0.017470</td>\n",
       "      <td>0.010746</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.027492</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>张家口市</td>\n",
       "      <td>3.680067</td>\n",
       "      <td>1.204148</td>\n",
       "      <td>2.334679</td>\n",
       "      <td>1.417932</td>\n",
       "      <td>1.110539</td>\n",
       "      <td>0.612830</td>\n",
       "      <td>0.439005</td>\n",
       "      <td>2.174117</td>\n",
       "      <td>inf</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.005701</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>承德市</td>\n",
       "      <td>2.454110</td>\n",
       "      <td>0.976658</td>\n",
       "      <td>1.627907</td>\n",
       "      <td>6.065174</td>\n",
       "      <td>3.048455</td>\n",
       "      <td>0.362100</td>\n",
       "      <td>0.222737</td>\n",
       "      <td>1.148757</td>\n",
       "      <td>2.144588</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.004509</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10 rows × 353 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  Cities       北京市        天津市       石家庄市        唐山市       秦皇岛市        邯郸市  \\\n",
       "0    北京市       inf  11.855277   5.725960  12.749903   9.573570   6.598759   \n",
       "1    天津市  4.458202        inf   1.806775  20.849449   5.471267   5.470594   \n",
       "2   石家庄市  2.230990   1.815094        inf   3.832100   3.980366  11.405949   \n",
       "3    唐山市  1.755526   6.625320   2.427825        inf  23.561492   1.150411   \n",
       "4   秦皇岛市  0.743911   1.104723   0.999676  12.371138        inf   0.347803   \n",
       "5    邯郸市  3.058310   3.685172   9.266465   2.033850   1.596207        inf   \n",
       "6    邢台市  1.784559   1.592507  17.716423   1.595459   1.481331  17.228512   \n",
       "7    保定市  7.296308   3.120498  14.208043   3.131967   2.768861   2.910327   \n",
       "8   张家口市  3.680067   1.204148   2.334679   1.417932   1.110539   0.612830   \n",
       "9    承德市  2.454110   0.976658   1.627907   6.065174   3.048455   0.362100   \n",
       "\n",
       "         邢台市        保定市       张家口市  ...      和田地区  伊犁哈萨克自治州      塔城地区  \\\n",
       "0   4.861948  27.363624  37.854040  ...  0.729644  0.424284  0.146950   \n",
       "1   1.854504   3.523714   3.469709  ...  0.527886  0.085204  0.012593   \n",
       "2  30.187961  15.307815   5.608304  ...  0.000000  0.081481  0.003760   \n",
       "3   0.838385   1.560130   2.586040  ...  0.000000  0.000000  0.000000   \n",
       "4   0.232399   0.585399   0.724584  ...  0.000000  0.000000  0.000000   \n",
       "5  18.648848   2.832027   1.373215  ...  0.000000  0.010470  0.000000   \n",
       "6        inf   2.246571   1.065858  ...  0.009354  0.031075  0.004738   \n",
       "7   2.589855        inf   5.916331  ...  0.000000  0.007716  0.005211   \n",
       "8   0.439005   2.174117        inf  ...  0.000000  0.000000  0.000000   \n",
       "9   0.222737   1.148757   2.144588  ...  0.000000  0.000000  0.000000   \n",
       "\n",
       "      阿勒泰地区      石河子市      阿拉尔市      五家渠市   香港特别行政区   澳门特别行政区  Unnamed: 352  \n",
       "0  0.244898  0.284529  0.192485  0.109344  1.716803  1.431816           NaN  \n",
       "1  0.069344  0.005139  0.027000  0.026867  0.228988  0.325382           NaN  \n",
       "2  0.000000  0.075537  0.006981  0.044484  0.012442  0.007495           NaN  \n",
       "3  0.000000  0.012455  0.000000  0.000000  0.000000  0.074082           NaN  \n",
       "4  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000           NaN  \n",
       "5  0.000000  0.040079  0.000000  0.000000  0.000000  0.000000           NaN  \n",
       "6  0.000000  0.034814  0.009440  0.009669  0.000000  0.000000           NaN  \n",
       "7  0.017470  0.010746  0.000000  0.027492  0.000000  0.000000           NaN  \n",
       "8  0.000000  0.000000  0.000000  0.005701  0.000000  0.000000           NaN  \n",
       "9  0.000000  0.000000  0.004509  0.000000  0.000000  0.000000           NaN  \n",
       "\n",
       "[10 rows x 353 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df=pd.read_csv('Pij_BAIDU.csv',encoding='gbk')\n",
    "df.head(10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00000000e+00, 4.45820151e-02, 2.23099039e-02, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [1.18552769e-01, 0.00000000e+00, 1.81509430e-02, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [5.72596015e-02, 1.80677508e-02, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       ...,\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [1.71680253e-02, 2.28988485e-03, 1.24418856e-04, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 2.05991973e-02],\n",
       "       [1.43181560e-02, 3.25381944e-03, 7.49467787e-05, ...,\n",
       "        0.00000000e+00, 1.74395725e-02, 0.00000000e+00]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cities = {d:i for i,d in enumerate(df['Cities'])}\n",
    "pij2 = np.zeros([len(nodes), len(nodes)])\n",
    "for k,ind in cities.items():\n",
    "    row = df[k]\n",
    "    for city,column in cities.items():\n",
    "        i_indx = nodes.get(city, '-1')\n",
    "        if i_indx < 0:\n",
    "            print(city)\n",
    "        j_indx = nodes.get(k, '-1')\n",
    "        if j_indx < 0:\n",
    "            print(k)\n",
    "        if i_indx >=0 and j_indx >= 0:\n",
    "            pij2[j_indx, i_indx] = row[column] / 100\n",
    "            if i_indx == j_indx:\n",
    "                pij2[i_indx, j_indx] = 0\n",
    "pij2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 将pij1和pij2合并\n",
    "合并思路，如果ij元素中pij2有值就用pij2的，否则将这一行所有的未被替换的pij1的值重新按照这些非零元素的数值大小计算占比，然后用1-$\\sum_{j}p_{ij2}$的总值乘以这个占比，从而得到新的数值，保证$\\sum_{j}p_{ij}=1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bools = pij2 <= 0\n",
    "pij = np.zeros([pij1.shape[0], pij1.shape[0]])\n",
    "for i in range(pij1.shape[0]):\n",
    "    row = pij1[i]\n",
    "    bool1 = bools[i]\n",
    "    values = row * bool1\n",
    "    if np.sum(values) > 0:\n",
    "        ratios = values / np.sum(values)\n",
    "        sum2 = np.sum(pij2[i, :])\n",
    "        pij[i,:] = (1 - sum2) * ratios + pij2[i, :]\n",
    "np.sum(pij,1) #验证一下是否归一化\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('pij.pkl', 'wb')\n",
    "pickle.dump(pij, f, pickle.HIGHEST_PROTOCOL)\n",
    "f.close()\n",
    "f = open('nodes.pkl', 'wb')\n",
    "pickle.dump(nodes, f, pickle.HIGHEST_PROTOCOL)\n",
    "f.close()\n",
    "\n",
    "f = open('city_info.pkl', 'wb')\n",
    "pickle.dump(city_properties, f, pickle.HIGHEST_PROTOCOL)\n",
    "f.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 如果没有原始数据，可以从这里开始\n",
    "加载数据:pij, nodes, city_properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('pij.pkl', 'rb')\n",
    "pij = pickle.load(f)\n",
    "f.close()f = open('nodes.pkl', 'rb')\n",
    "nodes = pickle.load(f)\n",
    "f.close()\n",
    "f = open('city_info.pkl', 'rb')\n",
    "city_properties = pickle.load(f)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 读取病例数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reader=pd.ExcelFile('city_cases.xlsx')\n",
    "df = reader.parse(\"sheet1\")\n",
    "df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_cases_cities = list(set(df['城市']))[1:]\n",
    "all_cases_cities = all_cases_cities + ['北京市','天津市','重庆市','上海市']\n",
    "#print(all_cases_cities)\n",
    "wuhan=df.loc[df['城市']=='武汉市',['新增确诊病例','确诊/出院','公开报道时间','新增治愈出院数','新增死亡数']]\n",
    "dates = list(wuhan['确诊/出院'])\n",
    "sorted_dates = np.sort(dates)\n",
    "first_date = sorted_dates[0]\n",
    "first_cases = 39\n",
    "print(first_date, first_cases)\n",
    "\n",
    "all_cases = {}\n",
    "\n",
    "for city in all_cases_cities:\n",
    "    subset = df.loc[df['城市']==city,['新增确诊病例','确诊/出院','公开报道时间','新增治愈出院数','新增死亡数','省份']]\n",
    "    zhixia = False\n",
    "    if len(subset)==0:\n",
    "        subset = df.loc[df['省份']==city[:-1],['新增确诊病例','确诊/出院','公开报道时间','新增治愈出院数','新增死亡数','省份']]\n",
    "        zhixia = True\n",
    "    new_cases = np.array(subset['新增确诊病例'])\n",
    "    cued_cases = np.array(subset['新增治愈出院数'])\n",
    "    die_cases = np.array(subset['新增死亡数'])\n",
    "    dates = list(subset['确诊/出院'])\n",
    "    dates1 = list(subset['公开报道时间'])\n",
    "    days = []\n",
    "    for i,dd in enumerate(dates):\n",
    "        if pd.isna(dd):\n",
    "            dd = dates1[i]\n",
    "        if not pd.isna(dd):\n",
    "            days.append(int((dd - first_date)/ np.timedelta64(1, 'D')))\n",
    "    sorted_days = np.sort(days)\n",
    "    indx = np.argsort(days)\n",
    "    infected = np.cumsum(new_cases[indx])\n",
    "    cued = np.cumsum(cued_cases[indx])\n",
    "    death = np.cumsum(die_cases[indx])\n",
    "    if len(sorted_days)>0:\n",
    "        if zhixia:\n",
    "            all_cases[list(subset['省份'])[0]+'市']= (sorted_days, infected, cued, death)\n",
    "        else:\n",
    "            all_cases[city] = (sorted_days, infected, cued, death)\n",
    "for case in all_cases.values():\n",
    "    xx = case[0]\n",
    "    yy = case[1]\n",
    "    if len(yy) > 0:\n",
    "        plt.semilogy(xx,yy,'o-')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wuhan\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_cases\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #读取工作簿和工作簿中的工作表\n",
    "# reader=pd.ExcelFile('Wuhan_nCoV202001241408.xlsx')\n",
    "# df = reader.parse(\"汇总表格\")\n",
    "# all_cases_cities = list(set(df['地级行政单位']))[1:]\n",
    "# wuhan=df.loc[df['地级行政单位']=='武汉',['病例数','发病日期','确诊日期','地级行政单位']]\n",
    "# dates = list(wuhan['发病日期'])\n",
    "# first_date = dates[0]\n",
    "# print(first_date)\n",
    "\n",
    "# all_cases = {}\n",
    "# for city in all_cases_cities:\n",
    "#     subset = df.loc[df['地级行政单位']==city,['病例数','发病日期','确诊日期']]\n",
    "#     a = np.array(subset['病例数'])\n",
    "#     dates = list(subset['发病日期'])\n",
    "#     dates1 = list(subset['确诊日期'])\n",
    "\n",
    "#     days = []\n",
    "#     for i,dd in enumerate(dates):\n",
    "#         #if pd.isna(dd):\n",
    "#         #    dd = dates1[i]\n",
    "#         if not pd.isna(dd):\n",
    "#             days.append(int((dd - first_date)/ np.timedelta64(1, 'D')))\n",
    "#     sorted_days = np.sort(days)\n",
    "#     indx = np.argsort(days)\n",
    "#     infected = np.cumsum(a[indx])\n",
    "#     if len(sorted_days)>0:\n",
    "#         all_cases[city] = (sorted_days, infected)\n",
    "# for case in all_cases.values():\n",
    "#     xx = case[0]\n",
    "#     yy = case[1]\n",
    "#     if len(yy) > 0:\n",
    "#         plt.semilogy(xx,yy,'o-')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 有效距离"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了计算有效距离，需要重新修改矩阵乘积的法则，基本思路就是在计算概率的时候，要把加法变为求大，计算距离的时候，要求小"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_matrix_production(matrix, matrix1):\n",
    "    #特殊的矩阵乘积，将加法转换为取大\n",
    "    sz = matrix.shape[0]\n",
    "    output = np.zeros([sz, sz])\n",
    "    for i in range(sz):\n",
    "        for j in range(sz):\n",
    "            row = matrix[i, :]\n",
    "            col = matrix1[:, j].reshape(-1)\n",
    "            prod = row * col\n",
    "            ele = np.amax(prod)\n",
    "            output[i,j]=ele\n",
    "    return output\n",
    "def eff_distance(prob):\n",
    "    #输入迁移概率矩阵，计算有效距离\n",
    "    sz = prob.shape[0]\n",
    "    prod = prob\n",
    "    distance = np.ones([sz, sz]) - np.log(prod)\n",
    "    for i in range(1, sz - 1):\n",
    "        flushPrint(i / sz)\n",
    "        prod = my_matrix_production(prod, prob)\n",
    "        dist = i + 1 - np.log(prod)\n",
    "        distance = np.minimum(distance, dist)\n",
    "    return distance\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 利用概率转移矩阵，计算有效距离"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "effective_distance = eff_distance(pij)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#outfile = open('effective_distance.npy', 'w')\n",
    "outfile = TemporaryFile()\n",
    "np.save(outfile, effective_distance)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 计算出武汉到每个城市的有效距离"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "距离排序并显示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eff_dist = copy.deepcopy(effective_distance)\n",
    "for i in range(len(nodes)):\n",
    "    eff_dist[i,i]=np.inf\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rowid = nodes['武汉市']\n",
    "sorted_distance = np.sort(eff_dist[rowid,:])\n",
    "sorted_cities = np.argsort(eff_dist[rowid,:])\n",
    "sorted_flux = np.sort(pij[rowid,:])[::-1]\n",
    "sorted_flux_cities = np.argsort(pij[rowid,:])[::-1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_in_nodes(ind, nodes):\n",
    "    key = -1\n",
    "    for k,v in nodes.items():\n",
    "        if v == ind:\n",
    "            key = k\n",
    "            break\n",
    "    return key\n",
    "firsts=[]\n",
    "firsts_dis=[]\n",
    "city_names_dis=[]\n",
    "for n, i in enumerate(sorted_cities):\n",
    "    city_name = find_in_nodes(i, nodes)\n",
    "    itm = all_cases.get(city_name,[])\n",
    "    if len(itm)==0:\n",
    "        itm = all_cases.get(city_name[:-1],[])\n",
    "    if len(itm)>0:\n",
    "        firsts.append(itm[0][0])\n",
    "        firsts_dis.append(sorted_distance[n])\n",
    "        city_names_dis.append(city_name)\n",
    "    print(city_name, sorted_distance[n])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n, i in enumerate(sorted_flux_cities):\n",
    "    city_name = find_in_nodes(i, nodes)\n",
    "    print(city_name, sorted_flux[n])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(eff_dist[rowid,:],pij[rowid,:],'o')\n",
    "plt.xlabel('Effective Distance')\n",
    "plt.ylabel('pij')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 与流量排序做比较"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,name in enumerate(city_names_dis):\n",
    "    print(i,name,firsts[i],firsts_dis[i])\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 有效距离和首感日期的比较"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "yy=np.array(firsts)*1.0\n",
    "xx=np.array(firsts_dis)\n",
    "bools=(xx!=inf)#&(xx>5.5)\n",
    "plt.figure(figsize=(15,10))\n",
    "\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(xx[bools],yy[bools])\n",
    "#slope1, intercept1, _,_,_ = stats.linregress(yy[bools],xx[bools])\n",
    "bools1 = (xx!=inf)&(xx>5.5)\n",
    "slope1, intercept1, r_value2,_,_ = stats.linregress(xx[bools1],yy[bools1])\n",
    "print(r_value,r_value2)\n",
    "\n",
    "bools2 = (sorted_distance!=inf)&(sorted_distance>=11)\n",
    "plt.plot(sorted_distance[bools2], slope1 * sorted_distance[bools2] + intercept1, 'ko',alpha=0.5)\n",
    "plt.plot(firsts_dis,firsts,'o')\n",
    "plt.plot(xx[bools1],yy[bools1],'o')\n",
    "\n",
    "plt.plot(xx[bools],slope*xx[bools]+intercept,'r-',label='{:2.2f}x+{:2.2f}'.format(slope,intercept))\n",
    "#plt.plot(xx[bools],(1/slope1)*xx[bools]-intercept1/slope1,'-',label=\n",
    "#         '颠倒x、y轴拟合：{:2.2f}x{:2.2f}'.format(1/slope1,-intercept1/slope1))\n",
    "plt.plot(xx[bools],slope1*xx[bools]+intercept1,'-',\n",
    "         label='去掉湖北省内的城市：{:2.2f}x+{:2.2f}'.format(slope1,intercept1))\n",
    "\n",
    "zhfont1 = matplotlib.font_manager.FontProperties(fname='/Users/zhangjiang/Library/Fonts/SimHei.ttf', size=10)\n",
    "show_cities = sorted_cities[bools2]\n",
    "for n,i in enumerate(show_cities):\n",
    "    xcor =sorted_distance[bools2][n]\n",
    "    ycor = slope1 * xcor + intercept1\n",
    "    #print(xcor, ycor)\n",
    "    city_name = find_in_nodes(i, nodes)\n",
    "    plt.text(xcor,ycor,city_name,fontproperties=zhfont1,alpha=0.5)\n",
    "for i in range(13):\n",
    "    plt.text(firsts_dis[i],firsts[i],city_names_dis[i],fontproperties=zhfont1)\n",
    "plt.text(firsts_dis[17],firsts[17],city_names_dis[17],fontproperties=zhfont1)\n",
    "plt.text(firsts_dis[19],firsts[19],city_names_dis[19],fontproperties=zhfont1)\n",
    "plt.text(firsts_dis[24],firsts[24],city_names_dis[24],fontproperties=zhfont1)\n",
    "plt.text(firsts_dis[25],firsts[25],city_names_dis[25],fontproperties=zhfont1)\n",
    "\n",
    "plt.xlabel('Effective Distance',fontsize=14)\n",
    "plt.ylabel('Reported Date',fontsize=14)\n",
    "#plt.ylim([15, 50])\n",
    "plt.legend(loc='upper left', shadow=True, numpoints = 1,fontsize=10,prop = zhfont1 )\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 传播模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "传播模型如下：\n",
    "\n",
    "\n",
    "$\\partial j_n/\\partial t = \\alpha s_n j_n - \\beta j_n + \\gamma \\sum_{m\\neq n}P_{mn}(j_m-j_n)$\n",
    "\n",
    "\n",
    "\n",
    "$\\partial s_n/\\partial t = -\\alpha s_n j_n + \\gamma \\sum_{m\\neq n}P_{mn}(s_m-s_n)$\n",
    "\n",
    "\n",
    "$r_n = 1 - s_n - j_n$\n",
    "\n",
    "其中$j_n$是n城市感染比例；$s_n$是n城市疑似人群比例；$\\alpha$为感染系数；$\\beta$为康复系数；$\\epsilon$为平均每城市人口数；$\\gamma$为平均流动人口占总人口的比例；$P_{mn}$为从n城市到m城市的n城市流动人口占比。\n",
    "\n",
    "在数值求解中，$\\sum_{m\\neq n}P_{mn}(j_m-j_n)$的计算用了矩阵和向量的乘法，并做了简化，即$\\sum_{m\\neq n}P_{mn}(j_m-j_n)=\\sum_m P_{mn}j_m-j_n$，写成矩阵方程为：$j\\cdot P - j$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "\n",
    "def protect_decay(t, t0, eta, rate_time):\n",
    "    epsilon = 0.001\n",
    "    r = 2 * np.log((1-epsilon)/epsilon) / rate_time\n",
    "    x0 = t0 + rate_time/2\n",
    "    decay = eta / (1 + np.exp(r * (t - x0))) + 1 - eta\n",
    "    return decay\n",
    "    \n",
    "def diff(sicol, t, alpha, beta, gamma, eta, rate_time, protect_day, pijt, intervention):\n",
    "    sz = sicol.shape[0] // 2\n",
    "    js = sicol[:sz]\n",
    "    ss = sicol[sz:]\n",
    "    jterm = js.dot(pijt) - js\n",
    "    sterm = ss.dot(pijt) - ss\n",
    "    if intervention:\n",
    "        cross_term = alpha * protect_decay(t, protect_day, eta, rate_time) * js * ss\n",
    "    else:\n",
    "        cross_term = alpha  * js * ss\n",
    "    delta_i = cross_term - beta * js + gamma * jterm\n",
    "    delta_s = - cross_term + gamma * sterm\n",
    "    output = np.r_[delta_i, delta_s]\n",
    "    return output\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 跑一下正式数据\n",
    "\n",
    "其中各个参数取值按照Science论文中的数值，可以考虑带代入我国参数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 计算流量比例数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timespan = np.linspace(0, 200, 1000)\n",
    "js0 = np.zeros(len(nodes))\n",
    "ss0 = np.ones(len(nodes))\n",
    "js0[nodes['武汉市']] = float(first_cases)/float(city_properties['武汉市']['pop'])#1e-4\n",
    "r0=3.7\n",
    "beta = 1/10 #一个病患从确诊到治愈（死亡）假设平均需要10天\n",
    "alpha = r0 * beta\n",
    "\n",
    "gamma = flowing_ratio\n",
    "print(gamma)\n",
    "\n",
    "eta = (1-1.0/r0)\n",
    "print(eta)\n",
    "rate_time = 30\n",
    "#武汉封城1月23日，数据起始日：1月10日\n",
    "result = odeint(diff, np.r_[js0,ss0], timespan, args = (alpha, beta, gamma, eta, \n",
    "                                                        rate_time, 13, np.transpose(pij),False))\n",
    "\n",
    "cities = ['武汉市','孝感市','黄冈市','深圳市','荆州市','襄阳市','上海市','北京市','黄石市','广州市']\n",
    "\n",
    "plt.figure(figsize = (15,8))\n",
    "plot_time_span = 700\n",
    "for i in range(result.shape[1]//2):\n",
    "    plt.plot(timespan[:plot_time_span],result[:plot_time_span,i],alpha=0.1)\n",
    "colors = plt.cm.jet(np.linspace(0,1,len(cities)))\n",
    "for n,i in enumerate(cities):\n",
    "    cityname = i\n",
    "    if cityname[-1]=='市':\n",
    "        cityname1 = cityname[:-1]\n",
    "    itm = all_cases.get(cityname1,[])\n",
    "    if len(itm)==0:\n",
    "        itm = all_cases.get(cityname,[])\n",
    "    if len(itm)>0:\n",
    "        idx = nodes.get(cityname,-1)\n",
    "        if idx>=0:  \n",
    "            plt.plot(timespan[:plot_time_span],result[:plot_time_span,idx],color=colors[n],label=cityname,linewidth=2)\n",
    "            infected_number = itm[1] - itm[2] - itm[3]\n",
    "            plt.plot(itm[0], infected_number / city_properties[cityname]['pop'], 'o',color=colors[n])\n",
    "\n",
    "zhfont1 = matplotlib.font_manager.FontProperties(fname='/Users/zhangjiang/Library/Fonts/SimHei.ttf', size=16)\n",
    "plt.legend(loc='upper left', shadow=True, numpoints = 1,fontsize=10,prop = zhfont1 )\n",
    "plt.xlabel('Day')\n",
    "plt.ylabel('Population Ratio')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 建模干预的情况"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "建模思路是将原方程改为：\n",
    "\n",
    "$\\partial j_n/\\partial t = \\alpha \\xi(t,t_0,t_m,\\eta) s_n j_n - \\beta j_n + \\gamma \\sum_{m\\neq n}P_{mn}(j_m-j_n)$\n",
    "\n",
    "\n",
    "\n",
    "$\\partial s_n/\\partial t = -\\alpha \\xi(t,t_0,t_m,\\eta) s_n j_n + \\gamma \\sum_{m\\neq n}P_{mn}(s_m-s_n)$\n",
    "\n",
    "\n",
    "$r_n = 1 - s_n - j_n$\n",
    "\n",
    "其中，\n",
    "\n",
    "$\\xi(t,t_0,t_m,\\eta) = \\frac{\\eta}{1 + \\exp(\\lambda (t - t_0 - t_m/2))} + 1 - \\eta$\n",
    "\n",
    "这里，\n",
    "\n",
    "$\\lambda = 2 \\frac{\\log(\\frac{1-\\epsilon}{\\epsilon})}{t_m}, \\epsilon = 0.01$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#有干预情况：\n",
    "\n",
    "timespan = np.linspace(0, 200, 1000)\n",
    "js0 = np.zeros(len(nodes))\n",
    "ss0 = np.ones(len(nodes))\n",
    "js0[nodes['武汉市']] = float(first_cases)/float(city_properties['武汉市']['pop'])#1e-4\n",
    "r0=3.7\n",
    "beta = 1/10 #一个病患从确诊到治愈（死亡）假设平均需要10天\n",
    "alpha = r0 * beta\n",
    "\n",
    "gamma = flowing_ratio\n",
    "print(gamma)\n",
    "\n",
    "eta = 1-4.0*beta/r0\n",
    "print(eta)\n",
    "rate_time = 90\n",
    "#武汉封城1月23日，数据起始日：1月10日\n",
    "result = odeint(diff, np.r_[js0,ss0], timespan, args = (alpha, beta, gamma, eta, \n",
    "                                                        rate_time, 13, np.transpose(pij),True))\n",
    "\n",
    "cities = ['武汉市','孝感市','黄冈市','深圳市','荆州市','襄阳市','上海市','北京市','黄石市','广州市']\n",
    "\n",
    "plt.figure(figsize = (15,8))\n",
    "plot_time_span = 1000\n",
    "for i in range(result.shape[1]//2):\n",
    "    for k,v in nodes.items():\n",
    "        if v == i:\n",
    "            cityname = k\n",
    "    plt.plot(timespan[:plot_time_span],result[:plot_time_span,i]*city_properties[cityname]['pop']\n",
    "             ,alpha=0.1)\n",
    "colors = plt.cm.jet(np.linspace(0,1,len(cities)))\n",
    "for n,i in enumerate(cities):\n",
    "    cityname = i\n",
    "    if cityname[-1]=='市':\n",
    "        cityname1 = cityname[:-1]\n",
    "    itm = all_cases.get(cityname1,[])\n",
    "    if len(itm)==0:\n",
    "        itm = all_cases.get(cityname,[])\n",
    "    if len(itm)>0:\n",
    "        idx = nodes.get(cityname,-1)\n",
    "        if idx>=0:  \n",
    "            plt.plot(timespan[:plot_time_span],result[:plot_time_span,idx]*city_properties[cityname]['pop']\n",
    "                     ,color=colors[n],label=cityname,linewidth=2)\n",
    "            infected_number = itm[1] - itm[2] - itm[3]\n",
    "            plt.plot(itm[0], infected_number, 'o',color=colors[n])\n",
    "\n",
    "zhfont1 = matplotlib.font_manager.FontProperties(fname='/Users/zhangjiang/Library/Fonts/SimHei.ttf', size=16)\n",
    "plt.legend(loc='upper right', shadow=True, numpoints = 1,fontsize=10,prop = zhfont1 )\n",
    "plt.xlabel('Day')\n",
    "plt.ylabel('Population Ratio')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
